SPATIOSPECTRAL CONCENTRATION ON A SPHERE* 
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Abstract. We pose and solve the analogue of Slepian's time-frequency concentration problem 
on the surface of the unit sphere to determine an orthogonal family of strictly bandlimited functions 
that are optimally concentrated within a closed region of the sphere, or, alternatively, of strictly 
spacelimited functions that are optimally concentrated within the spherical harmonic domain. Such 
a basis of simultaneously spatially and spectrally concentrated functions should be a useful data 
analysis and representation tool in a variety of geophysical and planetary applications, as well as 
in medical imaging, computer science, cosmology and numerical analysis. The spherical Slepian 
functions can be found either by solving an algebraic eigenvalue problem in the spectral domain or by 
solving a Fredholm integral equation in the spatial domain. The associated eigenvalues are a measure 
of the spatiospectral concentration. When the concentration region is an axisymmetric polar cap 
the spatiospectral projection operator commutes with a Sturm-Liouville operator; this enables the 
eigenfunctions to be computed extremely accurately and efficiently, even when their area-bandwidth 
product, or Shannon number, is large. In the asymptotic limit of a small concentration region and 
a large spherical harmonic bandwidth the spherical concentration problem approaches its planar 
equivalent, which exhibits self-similarity when the Shannon number is kept invariant. 

Key words, bandlimited function, commuting differential operator, concentration problem, 
eigenvalue problem, multitaper spectral analysis, spherical harmonic 
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1. Introduction. In a classic series of papers published in the 1960's, David 
Slepian and his colleagues solved a fundamental problem in communications engineer- 
ing, namely that of optimally concentrating a signal in both the time and frequency 
domains 0; 113 IE3 13 • The orthogonal family of data windows, or tapers, that 
arise in this context, and their multi-dimensional extensions l2nH3(il | have been used as 
the basis for the multitaper method of spectral analysis |46l: l6 ij , and for the analysis 
and representation of data in a wide range of physical, computational and biomedical 
disciplines (e.g., geodesy, seismology, optics, information theory, neurology and speech 
recognition). Time- frequency and time-scale concentration operators have been stud- 
ied in more general sett ing s and a variety of one- and multi-dimensional geometries 
by several authors [e.g., H [U [13 113 Hi EH 113 • 

In this paper we consider the simultaneous spatial and spectral concentration 
of a real-valued function of geographical position on the surface of the unit sphere. 
The spherical multitapers that we derive here should be useful in a number of geo- 
physical and planetary data analysis applications; this is our primary motivation to 
undertake the present study. In particular, we note that physical properties, such as 
the thickness or elastic strength of a planetary lithosphere, can be estimated from 
the cross-spectral properties of the surface topography and associated gravitational 
field |63j . Such data are most commonly available as bandlimited spherical harmonic 
coefficients, measured by artificial satellites or spacecraft. In many if not most ap- 
plications, planetary curvature prohibits the use of locally fiat approximations p^ . 
Thus, the determination of spatially localized estimates of planetary properties re- 
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quires spatiospcctral localization methods that go beyond those available in the plane 
[e.g.,|50|- Single spherical windows or tapers have been developed and apphcd in a 



number of recent studies [e.g., [la [13; |29| |53|; however, these are neither optimally 
concentrated, nor as reliable as an orthogonal family of multitapcrs in the extraction 
of robust localized statistical information from bandlimited spherical data [H^ . 

Following an initial and extremely insightful analysis by Griinbaum and his col- 
leagues Slepian's concentration problem on a sphere has, to our knowledge, been 
revisited quite rarely, b y w orkers interested in geodesy 0, magnetic resonance imag- 
ing of the human brain |47l | and planetary spectral analysis . Each of these studies 
treats a special case. In this paper, we pose and solve the spherical spatiospcctral 
concentration problem in its most general form, discuss a number of numerical imple- 
mentation methods, and analyze the flat-earth asymptotic limit, in which the spherical 
concentration problem approaches the corresponding problem on a plane. 

2. Slepian concentration problem. We begin with a brief review of the one- 
dimensional, continuous-continuous, time-frequency concentration problem. The re- 
sults are well known so we shall just articulate them without providing any derivations; 
our only objective is to provide a template for the spherical concentration problem, 
which we consider in the remainder of the paper. We use t and w to denote time and 
angular frequency, respectively, and adopt a normalization convention in which a real 
time-domain signal f{t) and its Fourier transform F(lo) are related by 



-1 /"OO pOC 

(2.1) fit) = 7r F{uj)e'^Uco, F{oj)= f{t)e 



' dt. 



The specific problem considered by Slepian [5y| is that of optimally concentrating a 
strictly bandlimited signal g{t), with a spectrum G{lj) that vanishes for frequencies 
|w| > W , into a time interval \t\ < T. No such bandlimited signal g{t) can be com- 
pletely concentrated within a finite interval by virtue of the Hcisenbcrg uncertainty 
principle 0; 113 ; the optimally concentrated signal is considered to be the one with 
the least energy outside of the interval: 



T 

g\t)dt 

(2.2) A = —j^ = maximum. 

g^{t)dt 

Bandlimited signals g{t) satisfying the variational problem H2.2|l have spectra G{uj) 
that satisfy the frequency-domain convolutional integral eigenvalue equation 

GiLj')du;' ^ AG(a>), \u^\<W. 

-Kyw - uj ) 

A closely related problem is that of concentrating the spectrum H{uj) of a strictly 
timelimited function h{t)^ that vanishes for times \t\ > T, into a spectral interval 
|w| < W . The appropriate measure of concentration in this case is 

w 

\H{uj)\'^duj 

(2.4) A — = maximum. 

\H{Lj)\^da; 
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Timclimitcd signals h(t) whose spectra satisfy the variational problem H2.4|l them- 
selves satisfy the time-domain eigenvalue equation 

(2.5) r ""Z^'Zp Ht')dt'^Xh{t), \t\<T. 

Both equations H2.3|l and 1)2. 5|l have the same eigenvalues 1 > Ai > A2 > • • • > 0, with 
associated time-domain eigentapers gi(t) , g2{t) , . . . and hi{t) , h2{t) , . . . which coin- 
cide within the interval \t\ < T, and associated eigenspectra Gi(a;), G2(w), . . . and 
Hi{u}), i?2(w), . . . which coincide within the interval |a;| < W. 

A change of both the independent and dependent variables transforms both \2.'6\i 
and 1)2. 5|l into the same dimensionless eigenvalue equation: 

(2.6) / ' sinTWix^ x') ^ ^^^^^^ ^ ^ 
J -I i:{x-x) 

Equation (|2.6|) shows that the eigenvalues Ai, A2, . . . and suitably scaled eigenfunctions 
^/'i(x), -02 (x), . . . depend only upon the time-bandwidth product TW . The sum of the 
eigenvalues is related to this product by 

(2.7) M^tK^'I^. 

Because of the characteristic step shape of the eigenvalue spectrum j^O; 113 j this so- 
called Shannon number |46l | is a good estimate of the number of significant eigenvalues, 
or, roughly speaking, the number of signals /(t) that can be simultaneously well 
concentrated into a finite time interval \t\ < T and a finite frequency interval \iu\ < W. 

The integral operator acting upon ip on the left side of equation 1)2. 6|l commutes 
with a second-order differential operator, 

(2.8) V^^{l^x^)^~T^W^x^, 

ax ax 

which arises in the separation of the three-dimensional scalar wave equation in prolate 
spheroidal coordinates Because of this, it is also possible to find the scaled 

eigenfunctions tpi{x), 'ip2ix), ... by solving the Sturm-Liouville equation 

(2.9) ^^i.x^)^ + {x-T'W'x')i'^0. N<1, 
ax ax 

where x 7^ A is the associated eigenvalue. 

The bandlimited prolate spheroidal eigentapers may be chosen to be orthonormal 
over the infinite time interval \t\ < 00 and orthogonal over the finite interval \t\ < T: 

/oo 
gag(}dt = 5ai3 and / QaQii dt ^ X^Sap- 
-00 J-T 

Almost all of the above results can be extended to the analogous spatiospectral con- 
centration problem for functions defined on the surface of the unit sphere. As we shall 
see, this two-dimensional problem is enriched by the arbitrary shape of the region of 
spatial concentration. 
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3. Preliminaries. The geometry of the unit sphere = {f : ||f|| = 1} is de- 
picted in Figure ITTl We denote the colatitude of a point f by < 6* < tt and the 
longitude by < (/> < 27r, so that f = {9, (f) represents a geographical position on the 
sphere. The geodesic angular distance between two points f and f' will be denoted 
by A. where 

(3.1) cos A = f • f ' = cos 6 cos 6' + sin 6 sin 9' cos{(j) — (f)'). 

We use R to denote a region of H. of area A, within which we seek to concentrate a 
bandlimited function of position f . The region may consist of a number of unconnected 
subregions, R ~ RiU R2 U ■ ■ ■, and it may have an irregularly shaped boundary, as 
shown. The region complementary to R will be denoted by Q — R. 





Fig. 3.1. Sketch illustrating the geometry of the spherical concentration problem. Lower 
right shows an axisymmetric polar cap of colatitudinal radius O, treated in Section\^ The 
area of the region of concentration, R = Ri U R2 D ■ ■ ■ , is denoted by A. 



3.1. Spherical harmonics. Since we restrict attention to real-valued functions, 
we use real surface spherical harmonics Yim{r) = Yim{0, 4>) defined by foUl^ 



(3.2) 

(3.3) 
(3.4) 



Y, 



Im \ 



V2X;],„|(0) COSTO0 if -/ < TO < 
Xlo{0) ifTO = 

\/2Xim{0) sinrruj) if Q <m<l, 



xue) = {-ly 



2nv 





'{i-my: 




_{i + my._ 




\ l+7n 







P/,„(cos6'), 
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The quantity < Z < cx) is known as the angular degree of the spherical harmonic, 
and —l<m<lis its angular order. The I oo asymptotic wavenumber associated 
with a harmonic of degree I is yWTT) w / + 1/2 d; m. The function PimifJ-) 
defined in H3.4|l is the associated Legendre function of integer degree I and order m. 
The spherical harmonics Yimir) are eigenfunctions of the Laplace-Beltrami operator, 

(3.5) ^dl + cote de + {sin 0)-^dl, 

with associated eigenvalues + 1). Our choice of the multiplicative constants in 
equations H3.2|l and H3.3|l orthonormalizes the harmonics on the unit sphere: 



(3.6) / Yi„iYi>,n' dil = Su>S„im'- 

The corresponding fixed-order orthogonality relations for Xim {&) and Pim {p-) are 

r 1 

(3.7a) / XijnXi>,nSmede = —Su>, 

Jo 27r 

(3.7b) £p,„.p„„,,._J_|±^,,„. 

The integral of a Legendre polynomial Pi (/i) = Pio {fi) over a cap cos 6 < < 1 is 



(3.8) 



/ Pi d^L = —— (cos 9) - Pi+i (cos 6)] , , 



where P_i(/x) = 1, and the product of Legendre fimctions at the same argument is 



Xira{e)Xi,^{e) = (-1)'" J2 



{272 +1) {21 + 1) {21' + 1) 
An 

n—\l — l' 



where the arrays of indices are Wigner 3-j symbols [13 |39| . Of the numerous three- 
term recursion relations involving the associated Legendre functions and their deriva- 
tives, we shall make use of two in this paper, namely 

(3.10a) {21 + l)/iP,„ = {l- m + l)P,+i,™ + {I + m)Pi^i,^, 

dp 

(3.10b) (1 - m')^ ^{1 + 1)m^/™ ~{l~m + l)P,+i,™. 

Finally, there are two relations involving sums of products of Legendre functions 
evaluated at different arguments that are useful in the discussion that follows. The 
first is the well-known spherical harmonic addition theorem jl2j | 

(3.11) YUmm{f') = Pi{r ■ f'), 



7n — — l 
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and the second is the Legendrc version of the ChristofFel-Darboux identity |59|; |6C 

L 



[^i-^i')Y,{2l + l) 



l—m 



{I -my. 

(/ + m)! 



PiMPuniii') 



(3.12) 



(L-m + 1)! 



PL+l,m{p)PLm{p') ^ PLm{p-)PL+l,m{p-')] ■ 



[L + m)! 

An apphcation of L'Hopital's rule in equation (|3.12|) covers the case when /i = /i'. 

3.2. Functions on the sphere. Let /(f) be a real-valued, square-integrable 
function on the unit sphere fi. Any such function can be expanded as a series of 
spherical harmonics: 



(3.13) 



OC / „ 
1 n ™ J ^ ^ 



dn. 



1=0 m=-l 



Equations (|3.13|l are the spherical analogue of the one-dimensional Fourier transform 
pair (|2.1|l . The finite character of the unit sphere quantizes the colatitudinal and 
longitudinal "frequencies" < I < oo and —l<m<l. We use a sans serif f to denote 
the ordered column vector of spherical harmonic coefRcients: 

( \ 

(3.14) f= . 

V : / 

The norm of a function /(f) in the spatial domain will be denoted by 



(3.15) 



1/11^2 



n 



and the norm of its spectral-domain equivalent f will be denoted by 

oo / 



(3.16) 



ifiiL = E E f> 



2 

Im ' 



1=0 m=- 



llfll 



Using this notation, Parseval's relation can be written in the form ||/||q 
The power spectral density or variance per spherical harmonic degree I and per unit 
area of a function /(f) is defined by 



(3.17) 



if?) 



21 



1 ' 

m— — L 



2 

Im ' 



We use (5(f , f ') for the Dirac delta function on the sphere, with the rcphcation property 



(3.18) 



Sir,r')fir')dn' ^fir). 



We can write 5{r,r') in any of the alternative forms 

S{r,r') = {sin0)-^S{e-0')S{cl) -(/)') 
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OO / 



1=0 m=-l 



(3-19) =f:f^)^'(-r'). 



=0 



The degree variance of a delta function is constant; its power spectral density is white: 

(3.20) (<5f) = — for all < Z < OO. 

All 

3.3. Bandlimited and spacelimited functions. We are interested in two sub- 
spaces of the space of all squarc-intcgrable functions on the unit sphere f2. We use 

(3.21) 5i = {.g: (.gf) = for L < / < oo} 
to denote the space of all bandlimited functions, 

L I 

(3.22) 9 = Y11 mmYim. 

1=0 m=-l 

with no power above a maximal spherical harmonic degree L, and we use 

(3.23) Sr^ {h:h = Q inVL- R] 

to denote the space of all spacelimited functions h{v) that are strictly contained within 
a region R. The space Sr is infinite-dimensional, but the dimension of the space of 
bandlimited functions is 

L 

(3.24) dim^L = ^(2^ + 1) = (L + 1)2, 

since the ordered column vector of spherical harmonic coefficients 

/ 500 



(3.25) g = 



\ 9LL 



associated with a function g{v) of the form (|3.22|) has {L -f 1)^ entries. Spatial and 
spectral seminorms analogous to (|3.15ll and H3.16|) are defined as 



(3.26) 11/11?,^= f .fdn^ i|f|li = E E fL- 

•^^ ;=o m=-l 

Within the space Sr of spacelimited functions, is a norm, and within the space 

Sl of bandlimited functions, ||g|||, is a norm; more generally, however, both ||/|||. and 

2 
L 
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4. Concentration within an arbitrarily shaped region. The uncertainty 
principle 0; stipulates that no function can be both strictly spacclimited and 
strictly bandlimitcd, i.e., no /(f) can lie in both subspaces Sr and Sl simultaneously. 
The objective of this paper is to determine those bandlimitcd functions g(r) G Sl that 
are as well contained within a spatial region R as possible, and those spacelimited 
functions h{r) € Sr whose power spectrum is as well concentrated within the interval 
< Z < L as possible. As in the time-frequency case, these two spatiospcctral 
concentration problems will be shown to be each other's duals. For brevity, in the 
discussion that follows, we shall frequently use the abbreviations 

oo / oo LI L 

(4.1) EE =E EE -E- 

/— m— — l Im l—O ni — — l Irn 

4.1. Spatial concentration of a bandlimited function. To maximize the 
spatial concentration of a bandlimited function g{r) G Sl within a region R, we 
maximize the ratio of the (semi)norms: 



11.911?^ 


[ g^ dn 
Jr 






Ml 


[ g^dfl 




Jn 



(4.2) A = ^ = ^ = maximum 



The two-dimensional variational problem H4.2|l is analogous to the one-dimensional 
problem (|2.2|l . Here, as there, the quantity < A < 1 is a measure of the spa- 
tial concentration. Upon inserting the representation H3.22|l of g{r) into H4.2|) . and 
interchanging the order of summation and integration, we can express A in the form 



A = 



/ ( E f'™^'™ E ^''™''^'''"' I ' 

9l'm' 

^ \ Im I'm' / 

E^'™ E ( / ^'™^''™' ) -9'' 

Im I'm' ' 



E ^'"^ E ^Irnyi'm' dfl 

Im I'm' ^ 



91' 



(4.3) 



E E 

7n,l'm' Ql'm' 

Im I'rn' 

L ■ 

E 



where in the final step we have used the orthonormality relation and defined 

the quadruply indexed quantity 



(4.4) Dim,i'm' = / Yi,„yi,,n' dQ. 



R 
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Upon introducing the {L + 1)^ x (L + 1)^ matrix 

-Doo.oo • • • Dqo^ll 



(4.5) D 



Dll.oo ' ' ' Dll.ll 



with elements Dim,i'm', where < I < L and —l^^'m < I, we can rewrite equa- 
tion (|4.2() as a classical matrix variational problem ,22] in the space Sl- 

(A a\ \ 9^°9 

(4.D) A = — — = maxmium. 

g^g 

Column vectors g that render the Rayleigh quotient A in equation (|4.6|l stationary are 
solutions of the (L+l)^x(i + l)^ algebraic eigenvalue problem 

(4.7) Dg = Ag. 

Equation H4.7|l is the discrete spherical analogue of the one-dimension spectral-domain 
equation (|2.3|) . The matrix D is real, symmetric and positive definite, 

(4.8) D"^ = D and g"^Dg > for all g 7^ 0, 

so the {L -\- 1)^ eigenvalues A and associated eigenvectors g are always real j^^. We 
order the eigenvalues Ai, A2, . . . , A(i+i)2 and eigenvectors g^^, gj, . . . , g(L+i)2 so that 

(4.9) 1 > Ai > A2--- > A(i+i)2 > 0. 

Each spectral-domain eigenvector g^.^, a = 1, 2, . . . , (L -I- 1)^ gives rise to an associated 
bandlimited spatial eigenfunction (f ) , defined by equation (|3.22|) . The first inequal- 
ity in (|4.9|) is strict because no bandlimited function can be completely confined within 
a region i?, and the last inequality is strict because of the positive-definite character 
of the matrix D. 

The symmetry = D also guarantees that the eigenvectors gi,g2, . • . ,g(L+i)2 
are mutually orthogonal jl^l • We choose them to be orthonormal: 

(4.10) gIg/3 = ^Q/3 and g^Dg^ = AaJ^^. 

The associated spatial eigenfunctions gi(f), (72(f), . . . , .g(L+i)2(f) are in that case both 
orthonormal over the whole sphere 17 and orthogonal over the region R: 

(4.11) / gagpdV, = Sai3 and / gagp dil ^ XaSap- 
Jn J R 

The two spatial-domain relations in equation H4.11|l are equivalent to the correspond- 
ing matrix spectral relations l|4.10|l . and they are analogous to the one-dimensional 
orthogonality relations (|2.1Q(I . The eigenfunction 51(f) associated with the largest 
eigenvalue Ai is the member of the space Sl of bandlimited functions that is most 
spatially concentrated within the the region i?, the eigenfunction 32 (f) is the next 
best concentrated function in Sl orthogonal to 51(f) over both 17 and R, and so on. 
Written out in full using index notation, the matrix eigenvalue equation 14.7|l is 

L 

(4.12) Dlni.l'm'gi'm' = Agjm- 

I'm' 
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Upon multiplying equation (|4.12|) by Yimi^) and summing over all < Z < i and 
—I < m < I, the right side yields Xg{r), and the left can be manipulated as follows: 

^ ^ DlmMm'gi'rn' j i1,ri (? ) = X] X! ( / ^/m O^i'm' (? ') d^l' j grm'>^m(f) 
Im \l'm' J Im I'm' ^"^^ ^ 

(4.13) = / i5(f,f')5(f')'^^^', 



JR 

where in the final step we have defined the bandlimited Dirac delta function 

47r 



(4.14) D{r,r')^Yl E >lm(f)>z„.(f') = E ( ^) ^'(^ ' f')- 



Z=0 m=-i /=0 



The above derivation demonstrates that the spectral-domain matrix eigenvalue equa- 
tion 14.7|l is equivalent to the spatial-domain integral eigenvalue equation 

(4.15) [ D{r,r')g{r')dn' ^ \g{r), f e O. 

Jr 

Equation H4.15|l is a homogeneous Fredholm integral equation of the second kind, 
with a symmetric, separable kernel [l^l^. Upon inserting the representations H3.22|l 
and H4.14|l into equation (|4.15(l . we recover the matrix equation (|4.7(l . so that the 
spectral-domain eigenvalue problem for g and the spatial-domain eigenvalue problem 
for a bandlimited g(r) G Sl arc completely equivalent. 

In summary, we can find an orthogonal family of bandlimited eigcnfunctions that 
are optimally concentrated within a region R on the unit sphere either by solving 
the {L + 1)^ X {L + 1)^ matrix eigenvalue problem (|4.7|l for the spectral-domain 
eigenvectors 91,92, • ■ • ,9(l+i)2j or by solving the Fredholm integral equation H4.15|l 
for the associated spatial-domain eigcnfunctions 32, • ■ • , 5(l+i)2. Either method 
determines the optimally concentrated eigcnfunctions at all points f e fi, i.e., both in 
the region R, where they are concentrated, and in the complementary region il — i?, 
where they exhibit inevitable leakage. 

4.2. Spectral concentration of a spacelimited function. Instead of seeking 
to concentrate a bandlimited function g{r) S Sl within a spatial region R, we may seek 
to concentrate a spacelimited function h{r) G Sr within a spectral interval < I < L. 
A suitable measure of concentration is then the spectral norm ratio, analogous to the 
one-dimensional ratio (|2.4(l : 

TC\ \ \ L 1=0 m = -l 

4. lb) A = -- — = = maxmium. 

00 I 



,2 

1=0 m=-l 



Upon inserting the representation of the spherical harmonic expansion coefficients, 
(4.17) /i/,„= / hYi^dn, 
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and interchanging the order of summation and integration, we can rewrite the ra- 
tio H4.16|l in the form 



A 



(4.18) 



V / h{r)Yir„{r)dn f h{r')Yir„{r')dn' 

Irn 

V / h[v)Yunmdn f h{r')Yirn{r') dn' 
^^/j(f) (^Yirair)YUr')^ hir')dn'dQ 
j^J^h{v) {^Yi,^{v)Yu^{v')^ hir')dn'dQ 

[ [ h{r)D{r,r')h{r')dndn' 
Jr Jr 



/i2(f) dn 



IR 

where in the final step we have made use of rephcation property (|3.18() of the delta 
function H3.19|l and the definition (|4.14() of the kernel D{r,r'). Functions h{r) G Sr 
that render the Rayleigh quotient (|4.18|) stationary are solutions of the Fredholm 
integral eigenvalue equation 



(4.19) [ D{r,r')h{r')dn' = Xh{r), r e R. 

Jr 



IR 

Equation H4.19|l is the spherical analogue of the one-dimensional time-domain eigen- 
value equation H2.5|l . In fact, this equation for /i(f ) G Sr is identical to equation (|4.15() 
for g{r) € Sl- The only difference is that equation H4.15|l is applicable on the entire 
sphere ft, whereas the domain of equation H4.19|) is limited to the region R, within 
which h{r) ^ 0. Evidently, the eigenfunctions h{v) that maximize the spectral norm 
ratio (|4.16|) are identical, within the region R. to the eigenfunctions gij) that maxi- 
mize the spatial norm ratio H4.2|l : 

g[v) if f e i? 
otherwise. 



(4.20) h{v) 



Every one of the [L + 1)^ bandlimited eigenfunctions G 5l gives rise to a space- 
limited eigenfunction ha G Sr^ defined by the restriction (|4.2U|) . The associated 
eigenvalues Aq are a measure of both the spatial concentration of the bandlimited 
eigenfunctions within the region R and the spectral concentration of the spacelimited 
eigenfunctions within the interval < Z < L. The fractional spatial energy 1 — Aq. 
leaked by an eigenfunction gaij) to the region — i? is identical to the fractional 
spectral energy leaked by /!„(?) into the higher degrees L < I < oo. If we had started 
with the variational prescription (|4.1()|) rather than H4.2|) . we could have obtained 
the integral equation (|4.15|l governing a bandlimited function g{r) £ Sl by simply 
extending the domain of solution of equation (|4.19|) to the whole sphere il. 

The spacelimited eigenfunctions defined by equation (|4.20() are orthogonal but 
not orthonormal over the whole sphere Q and over the region R: 

(4.21) / hahpdn^ / hahp dfl = XaSafj- 

Jn Jr 
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The spherical harmonic cocfBcients /i;,„, with < I < oo and ~l < m < I, are given 
in terms of those of gi„i, with < Z < L and — I < rn < I, hy 

L I' 

(4.22) him = ^ ^ DimMm'gi'm', 

l'=Q m' = -l' 

which reduces to /i;,„ = Ap;,„ for < I < L. In addition to the {L + 1)^ eigenfimc- 
tions with associated non-zero eigenvalues Ai, A2, . . . , ^(L+iy^ equation (|4.19|) has an 
infinite-dimensional null space of cigenfunctions with associated eigenvalue A = 0. 
Every function h{r) that vanishes in CI ~ R and has no energy whatsoever in the 
interval < Z < L is a member of this null space. 

4.3. Significant and insignificant eigenvalues. The sum of the eigenvalues 
of the matrix D defined in equation H4.5|l is 

(L + l)^ L I 

(4.23) ^c,=ti-0 = Y,Yl ^w™- 

a—l t—0 m— — l 

Upon substituting for the diagonal matrix elements Dim._im from equation 1)4. 4|l 
and making use of the spherical harmonic addition theorem l|3.11(l . this simplifies to 



L I „ 

1=0 •'R \m=-l ) 



(4.24) =(L + i)2^, 

where in the final step we have used the identity P/(l) = 1 to express the result in 
terms of the surface area, A = dfi, of the region of concentration R. 

We can alternatively obtain the result H4.23|l by starting with the spatial eigen- 
value equation H4.15|l . The kernel D{r,r') can be expressed in terms of the spatial 
eigenfunctions 51, ff2j • ■ • i g{L+i)^j which constitute a basis for 5^, in the form 

{L+lf 

(4.25) Div,r')^ ^ 5c(f)5c(f')- 

a=l 

The representation (|4.25|) is the spherical analogue of Mercer's theorem j23|63- Alge- 
braically, we can regard the relation H4.25|l as having been obtained from the original 
representation (|4.14|) by an orthogonal transformation from one basis, Yim, Q < I < L 
and — / < TO < /, to another, ga, a = 1,2, . . . , {L + 1)^. Setting f' = f in equa- 
tion H4.25)l . and integrating over the region R, we obtain 

(4.26) / Div,v)dn^ V / glir)dQ^ V A„. 
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Alternatively, setting f ' = f in equation H4.14|l and integrating over R yields 

(4.27) l D(r,r) dn = ±i^^^ J^ P.H) dn = (L + i f ^, 

Comparing equations H4.23|l - I|4.24|l and (|4.26|) - H4.27|I . we see that we can write the 
sum of the eigenvalues (|4.23|l in any of the equivalent ways 

r A 

(4.28) N= > Aa=trD= / D{r,r) dn = {L + 1)^ —. 

— ' /r 47r 

a=l •' ^ 

The quantity N is the spherical analogue of the Shannon number 1)2.7(1 in Slepian's 
one-dimensional concentration problem. Eigenfunctions ga{i) that are well concen- 
trated within the region R will have associated eigenvalues Aq that are near unity, 
whereas those that are poorly concentrated will have associated eigenvalues that 
are near zero. If, as in the one-dimensional problem, the spectrum of eigenvalues 
Ai, A2, . . . , A(i+i)2 has a narrow transition band from values near unity to values near 
zero, then the total number of significant (A^ ~ 1) eigenvalues will be well approx- 
imated by the (rounded) sum H4.23|l . For this reason, we expect to be a good 
estimate of the number of significant eigenvalues. Roughly speaking, the spherical 
Shannon number (|4.28() is the dimension of the space of two-dimensional functions 
/(f) that are both approximately limited in the spectral domain to spherical harmonic 
degrees < I < L and approximately limited in the spatial domain to an arbitrarily 
shaped region R of area A |^ • 

Rather than seeking a bandlimited function g(f) £ Sl that is optimally concen- 
trated within a spatial region R, we could have decided to seek one that is optimally 
excluded from R, i.e., one that is optimally concentrated within the complementary 
region Q — R. In that case we would have sought to minimize rather than maximize 
the Rayleigh quotient ((4.2() . In fact, all that we have found are the bandlimited func- 
tions g{r) G Sl that render the eigenvalue A stationary, so we have actually solved the 
containment problem 1(4. 2|l and the exclusion problem simultaneously. The optimally 
concentrated eigenfunctions and the optimally excluded eigenfunctions are identical, 
but with the ordering indices reversed, i.e., the bandlimited function that is most 
excluded from R and most concentrated within fl — R is (7(^+1)2 (f), with the smallest 
associated eigenvalue X(^L+iy- Whenever the area A of the region i? is a small fraction 
of the area of the sphere, A <C 47r, there will be many more well excluded eigenfunc- 
tions with insignificant (Aq ~ 0) eigenvalues, than well concentrated eigenfunctions 
with significant (Aq, « 1) eigenvalues, i.e., N <^ {L + 1)^. 

The sum of the squares of the {L + 1)^ bandlimited eigenfunctions (/^(f) is a 
constant, independent of position f on the sphere f2. This is another consequence of 
Mercer's theorem (|4.25ll and the definition H4.14|) : 

(4.29) E'.9^(r)-^(M) = ^ = ^. 

a=l 

If the first N eigenfunctions gi, 172, ■ • ■ , 57V have eigenvalues near unity and lie mostly 
within i?, and the remainder g^+iT gN+2, ■ ■ ■ ,g(L+i)^ have eigenvalues near zero and 
lie mostly in 17 — _R, then we expect the eigenvalue-weighted sum of squares to be 



(4,30) Y. ^^'^(^) - E - „"'^ l^^^ 

a—1 a—1 ^ 
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The terms with N+l<a<{L + 1)^ should be neghgible, so it is immaterial 
whether they are included in the sum H4.30|l or not. Taken together, the first N 
orthogonal eigenfunctions ga-,a = 1,2, ...,7V, with significant eigenvalues Aq. ~ 1, 
provide an essentially uniform coverage of the region R. This is really the essence of 
the spatiospectral concentration problem; the number of degrees of freedom is reduced 
from dim^L = [L + 1)^ to the Shannon number N = {L + 1)'^A/{Att). 

4.4. Abstract operator formulation. We conclude this section on the con- 
centration problem for an arbitrarily shaped region by reiterating the above results 
using an abstract operator notation. We use H to denote the operator that acts 
upon square-integrable functions /(f) in the spatial domain to produce the associ- 
ated infinite-dimensional column vectors f of spherical harmonic coefficients fim in 
the spectral domain, and we use 7i~^ to denote its inverse, so that 

(4.31) Tif = f and Ti'^i = f. 

We also introduce two projection operators, TZ and C, which project onto the space 
Spj of spacclimited functions and the space of bandlimited functions, respectively. 
The first of these acts to spatially restrict functions /(f) in the spatial domain. 



(4.32) 7^/(f) 



/(f) if f e i? 

otherwise, 



whereas the second acts to bandlimit column vectors in the spectral domain, 

(/oo \ / /oo 
/ r / 
Joooo / \ JLL 

Finally, we introduce a notation for the standard inner product in both domains: 
(4.34) {f,f')n= f ff'dn and (f, f )oo = f^f 

Parseval's relation can be written using this notation in the form (/, /')o = (f,f')oo- 
The spatial and spectral norms introduced in equations (|3.15|) and H3.16|l are given 
by 11/11^ = (/,/)a and ||f||2^ = (f,f)^. 

The spatial-concentration variational problem H4.2|l and the spectral-concentration 
variational problem (|4.16() can be written using this operator notation in the form 

(4.35a) A = — ■ ^ ^ — maxmium, 

The associated spectral and spatial-domain eigenvalue equations are 

(4.36a) {cnnn-^c){a) ^ 

(4.36b) (nn'^cmz)(nf) = x(nf), 

where we have made use of the facts that H and H^^ are each others' transposes, 
that both TZ and C are their own transposes, and that 7?.^ = TZ and ^ C. Equa- 
tions (|4.36|) are the operator equivalents of the algebraic eigenvalue equation H4.7|l and 
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the integral eigenvalue equation (|4.19|l . Any solution of equation H4.36a|l is a band- 
limited column vector of the form g = £ f, whereas any solution of equation (|4.36b|l 
is a spacelimitcd function of the form h = TZf. Both the spectral-domain operator 
>C7i7?.7i"^£ and the spatial-domain operator TifH^^CT-CR, are symmetric by inspection. 
Application of the operator product TC~^CTi acts to bandlimit an arbitrary function 
/, an operation referred to as spherical or uniform- resolution filtering, or as triangular 
truncation, in numerical analysis [3 l25t IS l5^. 

5. Concentration within an axisymmetric polar cap. We turn our atten- 
tion next to the special but important case in which the region of concentration is 
a circularly symmetric cap of colatitudinal radius Q, centered on the north pole, as 
illustrated in the lower right of Figure ITTI 



(5.1) 



{e ■.Q<9 <e}. 



In practical applications, the eigenfunctions that are optimally concentrated within 
such a polar cap can be rotated to an arbitrarily positioned circular cap on the unit 
sphere using standard spherical harmonic rotation formulae 0; 0; ^2 • 

5.1. Decomposition of the spectral-domain eigenvalue problem. The 

matrix elements (14.411 reduce in the case H5.1f) to 



(5.2) 



hn.l''. 



271 (5^7 



f 

Jo 



XimXi'm' smOdO. 



The Kronecker delta 5mm' renders the matrix D in equation 1)4. 5|1 block-diagonal: 



/ Do 



(5.3) 



-'±1 



'±1 



\ 



D±L / 



where every (L — m + 1) x (L — m + 1) submatrix D-i-,„ ^ Dq occurs twice as a 
result of the doublet degeneracy associated with ±m. Rather than solving the full 
{L+lY X {L+iy^ eigenvalue equation (|4.7)1 . we may solve a series of smaller eigenvalue 
equations, ^±mQ±m — ^±mQ±m-< o^"^ foi' each order ±m. 

We shall henceforth drop the identifying subscript and rewrite the fixed-order, 
spectral-domain algebraic eigenvalue problem D-|-,„g_|_^ = X±mQ±ra s-s, simply. 



(5.4) 

The (L - TO + 1) X (L - m + 
vector g G 5i are of the form 



Dg = Ag. 

1) matrix D and the [L 



l)-dimensional column 



(Dram ' ' ' D^aL 

\ \ 

Dhm ■ ■ ■ DlL 

where, for a particular order m. 



and g = 



V 9L 



(5.6) 



Dw = 27r / XimXvmSmOdd. 

Jo 
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The integral in equation H5.6|) can be evaluated with the aid of equations H3.8|l - H3.9|l : 



^ {21 + l){2l' + \) f I n l' \ f I n I 



Du- = (-1)™ VAfl^^Vfi^ ^ 



V 



/ V m -m 



(5.7) X [P„_i(cose) - P„+i(cose)]. 

We rank order the distinct L — m+l eigenvalues Ai, A2, . . . , Xt-m+i obtained by 
solving the fixed-order eigenvalue problem (|5.4|) so that 

(5.8) 1 > Ai > A2 > • • • > Al^,„+i > 0, 

and orthonormalize the {L — m + l)-diniensional eigenvectors g^, 92, ■ 

(5.9) g^g^ = Sap and g^Dg^ = XaS^p. 




0' 40* ISO'O' 40* IKrO* 40' laO'OC 40* 160* 



Fig. 5.1. Colatitudinal dependence of the first four spatial- domain eigenfunctions 
gi{0),g2{6),gz{6),g4,{9) of fixed order m = (top) to m = 4 (bottom). The radius of the 
polar cap is Q = 40°, and the maximal spherical harmonic degree is L = 18. Black curves 
show the concentration within the cap < 6 < 40°; grey curves highlight the leakage into the 
rest of the sphere, 40° <0 < 180°. The eigenvalues Ai,A2,A3,A4 expressing the goodness of 
the spatial concentration are indicated. The corresponding leakage in the spectral domain is 
illustrated in Fiaure ^ET^ 
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The associated bandlimited eigenfunctions <?i(^), 52(^)1 ■ ■ ■ ,5L-m+i(^^), defined by 



L 



(5.10) g{e) = Y,9iXim{e), 

l—m 

then satisfy the colatitudinal orthogonahty relations 

(5.11) 2tt / gagp sinO d9 ^ Sap and 27r / gagp smd dd ~ Xa6, 

Jo Jo 



a/3- 



The optimally concentrated spatial eigenfunctions g{r) for a given m are expressed in 
terms of the fixed-order colatitudinal eigenfunctions (|5.10|) by 

\/2g{9) cosni<j) if — L < m < 

(5.12) 5(^,0) = <( 9{0) if m = 

V2g{9) sin 7710 if < m < L. 

The four most optimally concentrated eigenfunctions gi{9), g2{0), g3{9), g4{6), for 
orders < m < 4 are plotted in Fieure lKTI The radius of the polar cap in this example 
is 6 = 40°, and the maximal spherical harmonic degree is L ~ 18. The first zonal 
(m ~ 0) cigcnfunction, gi{9), has no nodes within the cap 0° < 8 < 40°; the second, 
(72 (^), has one node, and so on. The non-zonal (m > 0) eigenfunctions all vanish at 
the north pole, = 0°. The first four zonal eigenfunctions, the first three m — 1 
and m = 2 eigenfunctions, and the first two m = 3 and m = 4 eigenfunctions are all 
very well concentrated (A > 0.9), whereas the fourth to = 3 and to = 4 eigenfunctions 
exhibit significant leakage (A < 0.1). The numerical methods used to compute the 
results shown here and in ensuing figures arc summarized in Appendix^ 

5.2. Decomposition of the spatial-domain eigenvalue problem. The inte- 
gral eigenvalue problem (|4.15|) in the spatial domain likewise decomposes into a series 
of fixed-order, one-dimensional Fredholm eigenvalue equations, 

(5.13) / D{9,9')g{9')sm9' d9' = \g{9), 0<9<tt, 







each with an TO-depcndent, separable, symmetric kernel, 

L 

(5.14) D{9, 9') = 27r ^ X,„(0)X,,„( 

l—m 



The results (|5.13|l and H5.14|l can either be obtained by multiplying the index form 
Dii'gv = Ag; of equation H5.4|l by Xi„i{9) and summing over m < Z < i, or by 
substituting the representation H5.10() - (|5.12() into equation H4.15|l . and using the or- 
thogonality of the longitudinal functions . . . , \/2cosm(/), . . . , 1, . . . , V^siuTOc/), . . . over 
the interval < (/> < 2ti. The (L — to+1) x (L — m + 1) matrix eigenvalue problem H5.4|l 
can in turn be derived starting from the separable Fredholm equation H5.13|l . so that 
the fixed-order spectral and spatial eigenvalue problems are completely equivalent. 
The fixed-order, spacclimited eigenfunctions 



(5.15) h{9) = 



g{9) if < 61 < e 
otherwise 
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satisfy an equation identical to H5.13|l . but onfy within the polar cap itself: 



(5.16) 



D{e,e')h{0')sme' de' = xh{9), o < e* < e. 



The eigenvalue A is a measure of both the spatial concentration of g{9) G Sl within 
< < Q and the spectral concentration of h{9) G Sr within < Z < L. The 
substitution /i = cos 9 converts equations (|5.13|l and ()5.16|l into 



(5.17a) 
(5.17b) 



cos 

1 



Di^i,fi')gi^i')df/ ^\g{n), -l<fi<l, 



D{n, fi') h{n') dfi' = Xh{^i), cose < ^ < 1. 



COS G 



d3 



9 = * 



-.iO' 



. a.' 



A— 1 ffMsie 



I 




u 
E 



n 



□ 



4,tHaB9 





la 127 



5S 127 



1fi 127 




Fig. 5.2. Squared coefficients hf of the first four spacelimited eigenf unctions hi{9), h2{6), 
h3{d), hiiO) of fixed order m = (top) to m — 4 (bottom). The radius of the polar cap is 
Q = 40°, and the maximal spherical harmonic degree is L — 18. Black curves show the power 
within the interval of concentration < I < 18; grey curves show the leaked power within 
19 < I < 127. Values of hf are in dB, normalized to zero at the individual maxima. The 
eigenvalues Ai,A2,A3,A4 specifying the goodness of the spectral concentration are indic ated. 
The corresponding bandlimited, spatial- domain eigenfunctions are shown in Fiaure l^~ll 
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The kernel D{fj,,fj,') can be simplified using the ChristofTel-Darboux identity H3.12|l : 

PL + ljn{lJ')PLm(lJ'') - PLmjiAPL + l.Tnii^') 

I 1 

fl- H 

where L'Hopital's rule eovers the case fi = /i'. The squared spherical harmonic coeffi- 
cients hf of the four best concentrated spacelimited eigenfunctions hi{6), h2{0), h^{d), 
hi{6) for < m < 4 are plotted versus I in Figure lS^ The cap radius 8 = 40°, band- 
width L = 18, and layout are the same as in Figure The maximal contribution 
to the Qfth zonal (m = 0) eigenfunction comes from the harmonic degree satisfying 
^Jl{r+Y) w Ti/{aQ)\ physically, this corresponds to fitting an integral number of 
asymptotic wavelengths ^-k j ^1(1 + 1) within the cap of diameter 20. 

5.3. Significant fixed-order eigenvalues. The number of significant eigenval- 
ues, or partial Shannon number, for each of the fixed-order eigenvalue problems (|5.4|l . 
(|5.13|l . (|5.16|l or H5.17|l can be computed using any of the equivalent formulae 

L — m+l L .0 „i 

(5.19) A^™= II^" = ^Ai=/ D{e,e)de= D{fi, ^i) d^i. 



(5.18) D{fi,^i') 



1)! 



2(L + m)! 




Fig. 5.3. Fixed-order eigenvalue spectra for an axisymmetric polar cap of radius O = 40°. 
The maximal spherical harmonic degree is L = 18. A different symbol is used to plot Xa versus 
rank a for each order < m < 5. The total number of fixed- order eigenvalues is L — m-\- 1; 
only the largest eight (\i through Ag ) are shown. Vertical gridlines and top labels specify the 
partial Shannon numbers Nm, rounded to the nearest integer. 



20 FREDERIK J. SIMONS, F. A. DAHLEN AND M. A. WIECZOREK 

We can write the final relation in equation H5.19|) using equation ()5.18|l in the form 

(5.20) TV™ = \ / [PUi,nPLrn - ^L^Pl+I,™] dfl, 

where the prime denotes differentiation with respect to fi. Figure 1^31 shows the fixed- 
order eigenvalue spectra for < to < 5. The cap radius is = 40° and the maximal 
spherical harmonic degree is L = 18, as in Figures fS.ll and 15.21 The partial Shan- 
non numbers N„i, computed by rounding equation 15.191) to the nearest integer, are 
shown. As in the case of the classical Slepian problem |30t [i^ Is^ the spectra have a 
characteristic step shape, showing significant (A « 1) and insignificant (A w 0) eigen- 
values separated by a narrow transition band. The partial Shannon number H5.19|l 
provides a good estimate of the number of well concentrated eigenfunctions; the first 
Nfn eigenfunctions all have a concentration factor exceeding A = 0.5. 

Once the L + 1 sequences of fixed-order eigenvalues (|5.8|l have been found, they 
can be resorted in accordance with the mixed-order ranking H4.9|) . The total number 
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Fig. 5.4. Reordered eigenvalue spectra (\a versus rank a) for axisymmetric polar caps of 
Q = 10°, 20°, 30°, 40° and a common maximal spherical harmonic degree L = 18. The total 
number of eigenvalues is (L -I- 1)'^ = 361; only Ai through Aeo are shown. Different symbols 
are used to plo t the various orders —11 < m < 11; the first six symbols are the same as those 
used in Fiaure lE~^ Vertical gridlines and top labels specify the rounded Shannon numbers. 
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of significant eigenvalues is 

L 

(5.21) N = Na + 2Y,Nm, 

m—l 

where the factor of two accounts for the ±m degeneracy. In Figure we show the 
reordered, mixed-m eigenvalue spectra for four different polar caps, with colatitudinal 
radii 9 = 10°, 20°, 30°, 40°; the maximal spherical harmonic degree is i = 18. The 
rounded Shannon numbers TV = 3, 11, 24, 42, lie in the middle of the steep, transitional 
part of the spectra, roughly separating the reasonably well concentrated eigensolutions 
(A > 0.5) from the more poorly concentrated ones (A < 0.5) in all four cases. 

We illustrate the pointwise sums of squares J2a 9ai^' 'P) ^^"^ J2a ^adai^i 4') i^i 
FigurcESl for polar caps of radu 9 = 10°, 20°, 30°, 40° and a bandwidth L ^ 18. The 
full unweighted sum of {L + 1)^ terms (dashed) is equal to N/A = {L + 1)^/ {4n) over 
the entire sphere 0° < 9 < 180°, in accordance with equation H4.29|l . The eigenvalue- 
weighted sums are, in contrast, concentrated within the polar cap < 6* < 9; solid 
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Fig. 5.5. Cumulative energy of the eigenfunctions concentrated within circularly sym- 
metric polar caps of colatitudinal radii B = 10° , 20° , 30° , 40° . The maximal spherical har- 
monic degree is L = 18; the Shannon numbers are N = 3,11,24,42. The sums of squares 
g'i{0, (f>) -\-gi{9, (/)) + ■■■ and \ig\{6, 4>) + A2(?2(^j 0) + ■ ' ■> '""c plotted versus colatitude along 
a fixed arbitrary meridian <j). Dashed lines show the full unweighted sums of (L + 1)^ terms, 
which converge to the constant value N/A over the entire sphere < < tt . Solid lines show 
the eigenvalue-weighted partial sums of N/2 and N terms, and the full sum o/(L + l)^ terms, 
which are concentrated within the polar cap < 9 < Q . 
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Fig. 5.6. Bandlimited eigenfunctions g[6,(f>) that are optimally concentrated within a 
circular cap of radius Q = 40°. Dashed circle denotes the cap boundary. The maximal 
spherical harmonic degree is L = 18, and the Shannon number is N = 42. Subscripts on 
the eigenvalues Xa specify the fixed-order rank. The eigenvalues have been resorted into a 
mixed-order ranking, with the best concentrated eigenfunction plotted on the top left and the 
32nd best on the lower right. Regions in which the absolute value is less than one hundredth 
of the maximal absolute value on the sphere are rendered in white; blue is positive, and red 
is negative. 
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lines of different shades of grey distinguish the sums carried up to the first Af/2, N, 
or aU {L + f )^ possible terms. Our expectation in equation H4.30|l is seen to be well 
realized. The uniformity of coverage over the region of concentration achieved using 
only the first N eigentapers is a key result, responsible for the success of the multitapcr 
method of spectral estimation |6J|. The utilization of a single, non-oscillatory data 
taper generally discards information near the boundary; however, this information is 
recovered by application of the other orthogonal tapers. The energy of the first N 
tapered data sets is proportional to the original energy of the data within the region 
of concentration, enabling spatially uniform extraction of statistical information while 
minimizing spectral leakage |6^]. Mark's change here. 

Finally, Figure shows a polar plot of the 32 eigenfunctions g{9, (j)), defined by 
equations (|5.12|) . that are optimally concentrated within a polar cap of radius Q = 40°. 
The maximal spherical harmonic degree is L = 18 and the Shannon number is iV = 42, 
as in Figures f5.1H5.3l The eigenvalue ranking is mixed-order, as in Figure and 
all degenerate v^cosmc/), V2sinTO(/) doublets are shown. The concentration factors 
1 < A < 0.849 and orders —5 < m < 5 of each eigenfunction are indicated. Blue and 
red colors represent positive and negative values, respectively; however, all signs could 
be reversed without violating the quadratic concentration criteria (|4.2|l and (|4.16(l . 

5.4. Commuting differential operator. The analysis of the one-dimensional 
time-frequency concentration problem was advanced considerably by Slepian's seren- 
dipitous discovery of the commuting prolate spheroidal differential operator (12.811 . Re- 
markably, there is also a differential operator, discovered by Griinbaum et al. [ijj , that 
commutes with the integral operator on the left side of equations H5.16|l and l|5.17b|l 
[see also . This second-order differential operator is given explicitly by 

(5.22) g = (cos e - cos 0)V^ + sine - L{L + 2) cos 6*, 

du 



where 



,2 d"^ — d 



(5.23) = — ^+cot6i— -m'(sin6l)-^ 

du du 

is the fixed-order Laplace-Beltrami operator H3.5(l . Rewritten in terms of /.t = cos 6*, 
the Griinbaum operator H5.22|l is 

771^ (cos Q — fl) 



(5.24) (? = A 



(cose-/i)(l-M');^ 
d^ 



-L{L + 2)^i- ^ 

1- fi 



Since commuting operators have the same eigenfunctions, wc can find the spacclim- 
ited, fixed-order eigenfunctions h(9) or /i(/i) by solving the differential eigenvalue 
equation gh = x/i, where x 7^ A is the associated Griinbaum eigenvalue. 

To confirm that g in equation (|5.24(l is the desired commuting differential opera- 
tor, we are required to show that 

(5.25) / gfj,D{fi, fi') h{fi') dfi' ^ / D{fi, fi')gf^'h{fi') du' , 

J cos J COS 

for an arbitrary spacelimited function h(n). There are two steps in the argument 
given by Griinbaum et al. |l9l |. The first is to show that the right side of equation 
H5.25|l can be rewritten as 

(5.26) [ D{fi,fi')gf,,h{fi')dfi' ^ j gf,,D{fi,fi')h{fi')dfi\ 

J COS 'J COS 
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and the second is to show that 

(5.27) g^Diti,fi') = g^,D{fi,ti'). 

The first result 15.26|l is a straightforward exercise in integration by parts; for any 
two functions and ?7(/x), it may be easily shown that 



C{Q'n)dii 



COS 



1 

COS 



(cos e - m)(1 - )C'?/ + L{L + 2)fi C?? 



(5.28) 



-m^(cos6 - /i)(l - /i^) ^Cri dfi 



cos 



where we have used a prime to denote d/dfi. Although it is not needed for the proof 
of equation H5.25|l . we note for future reference that the result H5.28|l is also valid if 
the integrations are carried out over the full interval — 1 < < 1. To verify the second 
result H5.27|l we make use of the relation V^Pim = + \)Pim to write 



1 ^ 

[Z(/ + l)-L(L + 2)](2; + l) 



l—m 



{l-m)\ 



{l + m)\ 



(5.29) 



i(l-/?)X](2Z + l) 

l—'rn 

i(i-M'2)x:(2/+i) 



{l-m)\ 



{I + m)\_ 

'{I - my. 



{i + my. 



-^PiMPimifJ-') 



Plm{^i)-J^Plm{^i')- 



An application of the Legendre derivative identity Ij3.10b|l transforms H5.29|l into 



\{^^-^^')j2[l^-{L + lY]{2l + l) 



l—m 
L 



{l-m)\ 



{l + m)\ 



PiMPimi^i') 



i(/.-/.')E(2/ + l) 



{l-m+\)\ 
{I + m)! 



(5.30) 



X [Pl+lm{^i)Plrn{^J■') - -Rjri (/^)-P/+l jn (m')] : 

and the ChristofFcl-Darboux identity (|3.12|) transforms equation H5.30|l into 



(5.31) 



1 — 771 

^ L I 

^{^l-^l')Y.{2l + l)Y,{2n+\ 



{l + m)\ 
[n — m)! 



PiMPimin') 



[n + m)! 



Pnm 



Interchanging the order of summation in the last line of equation H5.31|l wc obtain 
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+ J(M-M')E[^'-(i+in(2/ + i) 



2 



{I - my. 



{I + my. 

(5.32) - i(M - m') J2 + 1) \r^^] PuMPuMj^i^l + !)• 

-> ^ — ' [n + m)\ ^ — ' 



n—m 



The final sum over n is equal to (i + 1)^ — n^, so that the two terms in equation (|5.32|l 
cancel, (^^ — Qf^r)D{iJ,, fj!) = 0, and the commutation relation (|5.25() is confirmed. 

5.5. Griinbaum's equation. The above argument shows that we can compute 
the fixed-order, spacelimited, colatitudinal eigenfunctions hi{9), h2{0), . . . , /il_„i+i(6') 
either by solving the integral equation H5.13|l or by solving the differential equation 

(5.33) (cos e - cos 6')V2/i + sin 61 ^-L(L + 2) cos 61 /i = x/i, < 61 < 9. 

do 

The equivalent equation in terms of /i = cos 6 is in standard Sturm-Liouville form : 

(5.34) {ph'Y - qh + xph = 0, cos e < < 1, 

where the prime denotes differentiation with respect to fj., and where 

(5.35a) p(Ai) = (cose - /i)(l - ^i^), 

(5.35b) g(/i) = m'^ [1 - fj.^)-\iJ. - cos e) - L{L + 2) n, 

(5.35c) p{fi) = 1. 

Equation (|5.34|) must be solved subject to the requirement that h{^) remain finite at 
the endpoints /i = cos0 and fj, = 1. The associated variational problem is 

/ {ph'^ + qh^) d^i 

(5.36) X = 1 = minimum. 

/ ph^ dp. 

J cos 

All of the familiar Sturm-Liouville theorems apply. In particular, we know that equa- 
tion H5.34|l has a simple spectrum, with an infinite number of distinct eigenvalues 
Xi < X2 < ■ • - J having an accumulation point at infinity. The rank orderings of the 
eigenvalues xi: X2, ■ • ■ and the spatiospectral concentration factors Ai, A2, . . . , XL-m+i 
are reversed, so that the eigenfunction hi(9) associated with the numerically smallest 
eigenvalue xij which has no nodes in the polar cap < < O, is the most concentrated 
fixed-order eigenfunction; h2{0), which has exactly one node, is the next best concen- 
trated, and so on. Only the first L — m+1 eigenfunctions hi{9), h2{0), . . . , /iL_m+i(6') 
with non-zero eigenvalues Ai, A2, . . . , Al_„i+i are of interest in most applications. 
The remaining eigenfunctions /iL-m+2(^), ^L-m+3(^), • • • arc in the null space of the 
integral equation (|5.16|) . 



26 



FREDERIK J. SIMONS, F. A. DAHLEN AND M. A. WIECZOREK 



5.6. Commuting tridiagonal matrix. As in the case of (|5.13f) and H5.16|l . we 

are free to extend the domain of equation H5.33|l to the entire domain < 6 < n; in 
that case, the mrknown ftmction must be bandlimitcd rather than spacchmited: 

(5.37) (cose-cos6l)V2.9 + sin6'^-L(i + 2)cose'.g = x5, < 61 < tt. 

do 

Upon substituting the harmonic representation (|5.1U|) of g{9) into equation (|5.37|) . 
multiplying both sides by 2t: sin 6 Xiim{d), integrating over < ^ < tt, and invoking 
the orthogonality relation (|3.7b|) . we obtain the algebraic eigenvalue equation 

(5.38) Gg = xg, 
where 

/ G mm ^mL 

(5.39) G= : : 

\ Ghva • ■ ■ Gll 
is the (L — m + 1) x (L — m + 1) matrix with elements 

/•TT 

(5.40) Gu'=2tt XirnigXi,„^)sin9d0. 

Jo 

Equation (|5.38() is the spectral-domain version of the differential eigenvalue equa- 
tion 15.37|l just as equation (|5.4() is the spectral-domain equivalent of the integral 
eigenvalue equation (|5.13|) . 

As we have noted, the symmetry relation H5.28|l is valid even if the interval of 
integration is extended to < 6' < tt. as in equation H5.40|l . This shows that the 
Griinbaum matrix G is symmetric: 

(5.41) G = G"^. 

In addition, the matrices D and G commute with each other, 

(5.42) DG = GD, 

so they have identical eigenvectors. The index version of 



(5.43) ^'"G"'' = 2^ / Xhn{gXi,,n) sin0de = ^ GiMw 



The interior expression involving both the operator Q and integration over the re- 
gion of concentration < < © is the IV or VI element of the symmetric matrix 
product DG = (DG)^. Verification of the intermediate steps requires the use of the 
orthogonality relation Ij3.7bp . the operator identity (|5.27|) . and both the symmetry 
relation H5.28|l and its extension to the interval < 6* < tt. 

There are a number of ways to evaluate the matrix elements H5.40|l . perhaps the 
most straightforward of which is to make use of the relation V^X;™ = —1(1 -I- \)Xim 
and the Legendre identities H3.10|l . In fact, the Griinbaum matrix G is tridiagonal: 

(5.44a) Gh = -Z(/ + l)cose, 

(5.44b) Gu.. ^ G,,u = [1(1 + 2) - + 2)] ^j^^^y 
(5.44c) Gu> = otherwise, 

which is in agreement with the corresponding result given by Griinbaum et al. jl9l |. 




Fig. 5.7. Rank-ordered Griinbaum eigenvalues for polar caps of different colatitudinal 
radii, O = 10°, 20°, 30°, 40°, and a common maximal spherical harmonic degree L = 18. 
Separate sequences of eigenvalues xi, X2, ■ ■ ■ , XL-m+i for each angular order < m < L are 
connected by lines, with each sequence offset horizontally by its order m for clarity. 



The solution of equation H5.38|l offers a particularly attractive means of computing 
the eigenvectors Q G Sl and thus the optimally concentrated polar cap eigenfunctions 
g{0) G Sl, because it only requires the numerical diagonalization of a tridiagonal 
matrix G with analytically prescribed elements l|5.44|l . and a spectrum of eigenvalues x 
that is guaranteed to be simple. To illustrate this, we show the Griinbaum eigenvalue 
spectra for axisymmetric polar caps of radii O = 10°, 20°, 30°, 40° and a maximal 
spherical harmonic degree L = 18 in Figure 15.71 The ranked eigenvalues for every 
order < m < L are connected by lines, each sequence offset by its order to facilitate 
inspection. Thus, L + 1 eigenvalues xi, X2, ■ • ■ , Xl+i are plotted for m = 0, whereas 
a single eigenvalue xi is plotted for m = L. The roughly equant spacing and absence 
of a numerically troublesome tail of near-zero values, as in Figures 15.315. 41 is evident. 

Diagonalization of the Griinbaum matrix G enables the stable computation of 
optimally concentrated spherical Slepian functions that are bandlimited to very high 
angular degrees. To illustrate this, we show in Figure the first two zonal (m = 0) 
eigenfunctions gi{0), .92 (^) that are optimally concentrated within a polar cap of radius 
8 = 40°, for increasing bandwidths L = 10, 100, 300, 600. As the bandwidth increases, 
the eigenfunctions become increasingly concentrated near the pole = Q° . In any 
multitaper spectral application, a greater and greater fraction of any data near the 
boundary of the polar cap will be downweighted upon windowing; however, with 
increasing bandwidth and thus Shannon number, more and more well concentrated 
eigentapcrs become available, to enable the recovery of this lost information. 
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Fig. 5.8. Zonal (m = 0) eigenfunctions gi{0) (left) and §2(0) (right) that are optimally 
concentrated within an axisymmetric polar cap of colatitudinal radius O = 40°, for a range 
of maximal spherical harmonic degrees L — 10,100,300,600 (increasing shades of grey). 
Roundoff error prevents the accurate computation of the L = 300 and L — 600 eigenfunc- 
tions by double-precision numerical diagonalization of the matrix D; diagonalization of the 
Griinbaum matrix G overcomes this obstacle. 





Fig. 5.9. Optimally concentrated (top row) and optimally excluded (bottom row) zonal 
(m = 0) eigenfunctions, for a circular polar cap of colatitudinal radius O = 40° and a 
maximal spherical harmonic degree L = 18. The optimally excluded eigenfunctions cannot 
be accurately computed by double-precision diagonalization of the matrix D. Solution of the 
concentration problem for a polar cap of radius O — 140° gives rise to the same eigenfunctions 
gi{9),g2{0),g3{0),gi{9), ■ ■ ■ ,gi6{0),gi7{6),gi8{0),gi9{0), but in reverse order. 
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Concentration within a large rather than a small polar cap is another problem 
for which Griinbaum's method is ideally suited. To illustrate this, we show in Fig- 
ure 15.91 the first four (31,52,53, <?4) and the last four (.916, <7i7, 5i8, 519) zonal (m = 0) 
cigenf unctions, once again for a polar cap of radius O ~ 40° and a maximal spher- 
ical harmonic degree L — 18. As noted in Section [4.31 the eigenfunctions that are 
optimally excluded from the polar cap Q = 40° arc optimally concentrated within 
the much larger antipodal cap O 140°. The actual eigenvalues Aig, A17, Ais, A19 are 
many orders of magnitude smaller than the listed values, which simply represent the 
noise floor of our double-precision computations. The optimally excluded eigenfunc- 
tions can nevertheless be accurately computed by diagonalizing the tridiagonal matrix 
G; this is not otherwise possible in double precision arithmetic (see Appendix . 

5.7. Abstract operator formulation. The spatial-domain commutation rela- 
tion (|5.25(l can be expressed using the operator notation of Section ^31 as 

(5.45) {mr^cnn)g = g{mr^cnn). 

Since the Griinbaum operator Q acts only upon spacelimited colatitudinal functions 
h{9), it must satisfy 

(5.46) g = g'R = ng. 

Upon pre-multiplying equation H5.45|l by £7Y, post-multiplying it by Ti^^L, and mak- 
ing use of the relation H5.46(l and the fact that = C, we obtain 

(5.47) {cnnn^^c){cngn^^c) = {cngn'^c){cnnn'^c). 

Equation H5.47|l is the abstract operator formulation of the spectral-domain matrix 
commutation relation DG = GD. The operator equivalents of the spectral-domain 
eigenvalue equation H5.38() and the spatial-domain eigenvalue equation (|5.33() are 

(5.48a) {Cngn-^C)iC^) =xiCi), 

(5.48b) g(7^/) =x(7^/), 

where f{9) is an arbitrary colatitudinal function, neither bandlimited nor spacelim- 
ited. Because of the commutation relations (|5.45|) and (|5.47|) we are free to solve 
equations (|5.48|) rather than the fixed-order version of equations H4.36|l . to find the 
bandlimited eigenvectors g = £f and the spacelimited eigenfunctions h{9) ~ TZf{9). 

6. Continental Concentration. To illustrate the theory for an irregularly 
shaped region, we consider the spatiospectral concentration problem in six of the 
Earth's continental regions, listed in Table 16.11 in order of increasing size, together 
with their Shannon numbers for different bandwidtlis. The spectral analysis of data 
within either the Earth's continents or oceans has a number of applications in geodesy, 
geophysics and oceanography [e.g., |3 01; the spherical Slepian multitapers 
illustrated here should be ideally suited for this task. 

Figure 16.11 shows the eigenvalue spectra for five of the six regions (Greenland, 
Australia, North America, Africa and Eurasia) and four different bandwidths, L = 
6, 12, 18, 24, corresponding to (L + 1)^ = 49, 169, 361, 625 eige nfunctions each. The 
cutoff wavenumber associated with a bandwidth limit L is a/ L{L + 1) « L + 1/2 
divided by the Earth's radius 0; ll^; the cutoff wavelengths corresponding to the 
choices L ^ 6,12,18 and 24 are 6200,3200,2200 and 1600 kilometers, respectively. 
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Continental 


Area A/(47r) 




Shannon number N 




region 


in % 


L = 6 


L = 12 


L = 18 L 


= 24 


Greenland 


0.44 





1 


2 


3 


Australia 


1.50 


1 


3 


5 


9 


South America 


3.50 


2 


6 


13 


22 


North America 


3.98 


2 


7 


14 


25 


Africa 


5.78 


3 


10 


21 


36 


Eurasia 


9.98 


5 


17 


36 


62 



Table 6.1 

Areas, Shannon numbers, and bandwidths for the continental concentration problem. 



Only the largest of the Earth's continents, Eurasia, is sizable enough to exhibit at least 
one nearly perfectly concentrated eigenfunction for the smallest degree, L = 6, and the 
smallest region considered, Greenland, is too tiny to exhibit even a single eigenfunction 
with a concentration factor A near unity for the largest degree, L = 24. As in the case 
of a polar cap fFigure [5^ . the rounded Shannon numbers N = {L + 1)'^A/ (An) shown 
by the vertical dotted lines roughly separate the cigcnfunctions with concentration 
factors A > 0.5 from those with concentration factors A > 0.5. 
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Fig. 6.1. Eigenvalue spectra for five of the Earth's continental regions (Greenland, Aus- 
tralia, North America, Africa, Eurasia) and four different maximal spherical harmonic de- 
grees (L = 6, 12, 18, 24/ Vertical gridlines and five leftmost ordinate labels specify the rounded 
Shannon numbers N . Ordinates are truncated on the right; number to right of arrow is the 
total number of eigenvalues, (L + 1)^ = 49, 169, 361, 625. 
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Fig. 6.2. Bandlimited L = 18 eigenfunctions gi, g2, ■ ■ ■ , gi2 that are optimally concen- 
trated within the continent of Australia. The concentration factors Ai, A2, . . . , A12 are indi- 
cated; the Shannon number is N = 5. Order is left to right, top to bottom, as with English 
text. 



In Figurcs lF?^l6.3l and l6.4l wc show map views of the first twelve i = 18 eigenfunc- 
tions 51(f), g2(f), • ■ • 1512(f) that are optimally concentrated within Australia, North 
America and Africa. Blue colors denote positive values and red colors denote negative 
values (though, as we have noted, these could be reversed, since the sign of an eigen- 
function is arbitrary). Regions in which the absolute value is less than one hundredth 
of the maximal absolute value on the sphere are rendered in white. 

In the case of Australia (Figure l^^l the first five eigenfunctions are reasonably 
well concentrated within the continental boundaries (A5 = 0.607); however, the con- 
centration factors A diminish rapidly thereafter, so that gi2 is far more excluded than 
concentrated (A12 = 0.049). With a limiting bandwidth L = 18, and thus a cutoff 
wavelength of 2200 kilometers, it is only possible to concentrate = 5 orthogonal 
bandlimited eigenfunctions .91 , .92 , 53 , .94 , 55 into a continent which, across its north- 
south waist, is only about 1500 kilometers wide. 

This situation is much improved in the case of the North American continent 
fFigure l^T5|l . which has an area A that is 2.7 times larger than the area of Australia. In 
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Fig. 6.3. Bandlimited L = 18 eigenfunctions gi, g2, ■ ■ ■ , gi2 that are optimally concen- 
trated within the continent of North America. The concentration factors Xi, A2, . ■ . , A12 are 
indicated; the Shannon number is N = 14. Format is identical to that in Fiaure W^ 

fact, North America has = 14 reasonably well concentrated L = 18 eigenfunctions, 
of which only the first twelve are shown. The first eigenfunction, 171, is a roughly 
circular dome centered in the middle of the continent, as in the case of Australia. 
Subsequent orthogonal eigenfunctions (72,53, ■• ■ exhibit lobes in previously uncovered 
regions. The coastline of North America is more irregular than that of Australia; 
Quebec and the Northwest Territories of Canada are essentially uncovered until g^, 
western Alaska is not adequately covered until gg and gio, and Florida, Baja California 
and southern Mexico are only covered by gn and 512 at the expense of substantial 
leakage (An = 0.642, A12 — 0.596) outside of the continental boundaries. 

Africa (Figure [6.4(1 . which has an area A that is 3.9 times larger than that of 
Australia, has = 21 reasonably well concentrated L = 18 eigenfunctions, the twelfth 
of which has a concentration factor A12 = 0.887. Once again, gi is roughly circular, 
and the subsequent orthogonal eigenfunctions 52 , 53 , ■ • ■ successively cover previously 
uncovered regions. West Africa is uncovered by 171 and 172, but becomes reasonably 
well covered by (73 and 55; likewise, South Africa is uncovered until 174 and 55. Other 
geographical features become well covered by the increasingly oscillatory orthogonal 
eigenfunctions (e.g., Egypt by 177 and 512). 

Figure shows the eigenvalue- weighted sum of squares X^q, '^a5a(f ) of the band- 
limited L = 18 eigenfunctions of all six of the Earth's continents (excluding Antarc- 
tica). We find the eigenfunctions (71(f), g2(f), . . . , 5(^+1)2 (f) by diagonalization of the 
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Fig. 6.4. Bandlimited L = 18 eigenfunctions gi, g2, ■ ■ . , gi2 that are optimally concen- 
trated within the continent of Africa. The concentration factors Ai, A2, . . . , A12 are indicated; 
the Shannon number is N — 21. Format is identical to that in Fiaure W^ 

{L + 1)"^ X [L + 1)^ matrix (|4.5|l formed by smrmring the corresponding matrices 
Dsurasia + DAfrica + • • • of cach of the six Continents. The combined area of all six 
continents is A/{4:Tt) = 25.2%, and the Shannon number is N = 91; the partial sums 
of the first N/A, N/2 and A'^ terms, as well as the full sum of all (L + 1)^ — 361 terms, 
are shown. The ability of the first N eigenfunctions to provide uniform coverage of 
the target area is evident; as in Figure 1^31 the coverage is only marginally improved 
by adding the remaining, poorly concentrated (L + l)'^ — N = 250 terms. Due to their 
small size, Australia and Greenland do not appear until the 1 N/2 and 1 N 
partial sums, respectively. Even then, the coverage of Greenland is imperfect, an ex- 
pected consequence of the small Shannon number for Greenland (TV = 2 for a maximal 
spherical harmonic degree L ~ 18). 
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Fig. 6.5. Cumulative eigenvalue-weighted energy of the first N/A, N/2, N and all (L + l)'^ 
eigenfunctions that are optimally concentrated within the ensemble of continents Eurasia, 
Africa, North America, South America, Australia and Greenland. The maximal spherical 
harmonic degree is L = 18; the cumulative fractional area is A/{4:Tt) = 25.2%; the Sh annon 
number is N — 91. Darkest blue on color bar corresponds to the expected value of the 

sum, as shown. Regions in which the value is less than one hundredth of the maximal value 
on the sphere are rendered in white. 



7. Asymptotic scaling. As we have noted, the eigenvalues Ai, A2, . . . and suit- 
ably scaled eigenfunctions V'i(a;), '!/'(2;)2, ... of the original Slepian concentration prob- 
lem (|2.6|l depend only upon the Shannon number N = 2TW/tt. This Shannon-number 
scaling is the only important feature of the onc-dimcnsional problem that does not 
carry over to the spatiospectral concentration problem on a sphere. Fundamentally, 
this lack of scaling is a consequence of the fact that it is not possible to shrink or 
magnify a region, such as Africa, on a sphere O of fixed radius ||f |j = 1, while keeping 
the angular relationships among all of the interior points the same. Shannon-number 
scaling on a sphere is exhibited only asymptotically, in the limit 

A 

(7.1) y4->0, L^oo, with N=(L + lf— held fixed. 

An 

In that limit of a small concentration area A and a large bandwidth < I < L, the 
curvature of the sphere becomes negligible and the spherical concentration problem 
approaches the concentration problem on a plane. 

7.1. Hilb approximation and Poisson sum formula. Two results underlie 
the consideration of the flat-earth limit (|7.1() . which we undertake in this section. The 
first is Hilb's asymptotic approximation for the Lcgcndrc functions [5 Isl Elt l60| . 

(7.2) XUO)^{-iy''\f^^^^f^Jrn[{l + 1/2)0], 0<9«7r, 



SPATIOSPECTRAL CONCENTRATION ON A SPHERE 



35 



where J,„(a;) is the Bcssel function of the first kind, and the second is the truncated 
Poisson sum formula, 



L oo 

(7.3) ^/(Z + l/2)= ^ (-l)W f{k)e--'-^->^dk, 

1=0 s=-oo 

which is valid for an arbitrary continuous function f{x). To verify the relation H7.3|l 
we start with the Fourier series representation of /(x) on the interval < x < 27r, 

(7.4) ^^"^^=2^ E / fi^y-'^^~''^du. 

Letting x — > x + 2t:1 in equation (|7.4|) and summing over < I < L yields 

(7.5) E/(^ + 2^^) = ^ E / f{u)e^'^^--^du, 

1=0 s=-oo'^0 

where we have shifted the interval of integration for each term by 2t:1 and made use of 
the 27r periodicity of the exponential. Upon dividing the argument by 27r, substituting 
k — uH2'k), and setting x = tt, wc obtain the desired identity H7.3|l . 

7.2. Scaled integral equation for an arbitrary region. An application of 
both the Hilb approximation (|7.2|1 and the Poisson sum formula (|7.3() enables us to 
write the Fredholm kernel D{y, f') in equation H4.15|l in the form 

" ^21 



L+l 

2'Kisk 



1=0 



^ ^ X^(/ + l/2)Jo[(? + l/2)A] 



27r V sin A 

1=0 



(7.6) =^V^E(-1)V„ •^o(fcA)e— -fcdfc. 



s— — oo 



Upon substituting k = {L + l)p and taking the limit L oo, A — > 0, with the product 
LA held fixed, equation H7.6|l reduces to 

(7.7) D{A)^^—J^ Jo[{L + l)pA]pdp= — , 

where we have made the approximation A/sin A sa 1, and used the Riemann-Lebesgue 
lemma to eliminate the s 7^ terms involving the highly oscillatory factors 
g-27rjs(L+i)p^ In the limit x ^ the ratio Ji(x)/x ^ 1/2, so the A ^ limit of 
the kernel (|7.7|l is D{Q) = {L + l)^/(47r), guaranteeing that the Shannon number, or 
sum of the eigenvalues (|4.23|l . is still given in this asymptotic approximation by 

(7.8) ^ = / D{0)dn = (i + 1 

Jr 

To obtain a scaled version of equation (|4.15|l dependent only upon the Shan- 
non number N, we make use of the approximation (|7.7|) for the kernel _D(f, f'), and 
introduce the independent and dependent variable transformations 

(7.9) x=\/^f, x' = \/^f', V^(x)=5(f), ^(x')=ff(f'). 



2 A 

Att' 
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The scaled coordinates x, x' are the projections of the points f , f' G 51 onto a large 
sphere f2» of squared radius ||xp = An/ A. The geodesic distance between the scaled 
points X, x' G 17, and the differential surface area on are 

/ 47r Att 

(7.10) l|x-x'|| = W — A and dfl^ ^ — dQ. 

V -/i -/I 

Upon making the substitutions 17.9|1 - H7.1U|1 . equations (|4.15|1 and 17.7|1 reduce to 

(7.11) / Z?,(x,x')V(x)dr!: = A7A(x), 

where i?, is the projection of the region of concentration R onto the sphere Q^, and 

^7 19^ n( VA^^i(VlVllx-x'll) 

(7.12) D4x,x) = -— 



27r ||x-x'|| 

is the symmetric, A^-dependent Fredholm kernel. 

Equations (|7.11|I - H7.12|I are the spherical analogue of the one-dimensional scaled 
eigenvalue equation H2.6|l . The asymptotic eigenvalues Ai, A2, . . . and associated scaled 
eigenfunctions ?/'i(x), ?/;2(x), . . . depend upon the maximal degree L and the area A 
only through the Shannon number N = {L + l)^A/(47r). As in the case of equa- 
tions (|4.15(l and H4.19|l . we are free to solve H7.11|l - H7.12|l either on all of fl^, in which 
case the eigenfunctions ?/'i(x), ■02 (x), . . . are bandlimited, or only in the region of con- 
centration i?,, in which case they are spacelimited. It is readily verified that the 




4S* 3D- 135' 

angular distance A 



1£0» 0* 



43* 90* 133" 
angulaf distanced 



Fig. 7.1. Compa rison of the exact scaled kernels I7.14t with the flat-earth asymptotic 
approximation \7. j5| ) (black). The Shannon number h = 3, 10, 23, 40 is kept constant in 
each of the four panels, and the bandwidth used to compute the exact scaled kernels varies 
between L — 1 (worst fitting) and L = 100 (best fitting). 
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scaling has no effect upon the sum of the eigenvalues, inasmuch as 

(7.13) N = [ D,XO)dn,^^ I dn, = N. 

Jr., 47r Jji^ 

We show in Appendix that the scaled eigenvalue problem (|7.11|) - (|7.12|l is identical 
to that governing the two-dimensional concentration problem on a plane. 

The above considerations show that in the limit (|7.1|l . we expect the exact Fred- 
holm kernel H4.14|l . evaluated on 51* and normalized by its value at zero offset, 



(7.14) Di^A^A) 



^E(2^ + l)P.(cosyf A^, 



D{0) (L 

to be well approximated by the similarly normalized asymptotic kernel 

^^(A) _ 2Ji(V]VA) 



(7.15) 



D,{0) ViVA 



The quality of this asymptotic approximation to the kernel and the associated flat- 
earth scaling are illustrated in Figure ITTI In the four examples shown, with Shannon 
numbers iV = 3, 10, 23,40, the approximation is excellent even for angular distances 
as large as A w 135°, once the maximal spherical harmonic degree exceeds L = 3-4. 

7.3. Scaled eigenvalue equation for an axisymmetric polar cap. The flat- 
earth asymptotic version of the fixed-order colatitudinal eigenvalue problem 15.13|l can 
be obtained in two different ways; either by an application of the Hilb approxima- 
tion (|7.2|) and the Poisson sum formula H7.3|) to the kernel D(0,ff) given in equa- 
tion (|5.14|l . or by using the addition theorem for Bessel functions |27| . 

oo 

(7.16) Jo(fcA) = Jo{k9)Jo{ke') + 2 ^ J™(fc6l) J™(fc6i') cos m(</> - (j)'), 

m—l 

the representation H5.12|l of g{d,(f>), and the orthonormality of the longitudinal func- 
tions . . . , cosm0, . . . , 1, . . . , ^/2amm(f>, . . . over the interval < (/) < 27r to decom- 
pose equations (|4.15|) and H7.7|l into a series of individual eigenvalue problems, one 
for each order < m < L. Using either method, we find that equation H5.13|l can be 
approximated in the limit H7.1|l by 

(7.17) / Die,e')gie')e'de' ^Xg{e), 

where 







(7.18) D{e,e') = {L + if [ j^[{L + i)pe\j^[{L + i)pe']pdp. 

Jo 

It is convenient in the present instance to approximate the area of the small polar cap 
hy A = 27r(l — cos0) « ttQ^, and to introduce scaled coordinates that are slightly 
different from those in equations H7.9|l . namely 



(7.19) x = e/e, x' = e'/e, 4>{x)^g{e), ijj{x') = g{e'). 
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This leads to a scaled, fixed-order eigenvalue problem, 



(7.20) 



with an associated kernel 



D4x,x')i;{x')x'dx' = AV'(x), 



(7.21) 



D4x,x')=4:N / J„,{2VNpx) Jm{2^/Npx')pdp. 



whose eigenvalues Ai, A2, . . . and associated scaled eigenfunctions il)i(x), 'ip2{x)^ . . . de- 
pend upon the maximal spherical harmonic degree L and the cap radius O only 
through the small polar-cap Shannon number N ~ j{L + 1)^0^. 

Although the polar-cap scaling relations H7.20|l - H7.21|l are strictly valid only in the 
asymptotic limit i — > 00, ^ 0, the approximation is excellent even for moderate 
bandwidths L and sizable cap radii 8. For a fixed Shannon number N = 40, a maxi- 
mal degree in the range 25 < L < 40, and therefore a cap radius Q = /{L + 1) 
in the range 29° > 8 > 18°, the agreement between the fixed-order, scaled eigenfunc- 
tions is always within a few percent. In principle, the asymptotic results H7.20|l - H7.21|l 
would enable the determination of approximate polar cap eigenfunctions g(9) for vary- 
ing values of L and 8 by scaling a pre-computed catalogue of fixed- eigenfunctions. 
In practice, the construction and diagonalization of the tridiagonal Griinbaum ma- 
trix H5.44(l is so straightforward and efficient that it is preferable to simply compute 
the optimally concentrated eigenfunctions g{0) exactly. 

7.4. Asymptotic fixed-order Shannon number. The asymptotic approxi- 
mation to the number of significant eigenvalues associated with a given order m is 



N„ 



D<,{x, x) xdx 



'0 

4,N 



1 pi 



J^j (2 viV px^ pdpxdx 



2N 



Ji{2VN) + Jl^,{2VN) 



(7.22) 



(2m -f- 1) VA J„ (2 VA ) Jm+i (2 VA ) 

m 

l~J^{2VN)-2j2Jn{2VN) 



m 
2" 



?i=i 



The relationship H5.21|l between the total number A' of significant eigenvalues and the 
number A,„ associated with each order m is preserved in this asymptotic approxima- 
tion, inasmuch as 



(7.23) 



AT = Ao + 2 £ AT™ 

00 

j2(2VApq) + 2 JlpVNpq) 



rn—l 
1 pi 



'0 JO 

»1 pi 



m— 1 



p dp xdx 



AN I I pdpxdx = N. 
'0 Jo 
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In Figure 17.21 we compare the exact fixed-order Shannon numbers N„i , computed 
by Gauss-Legendre numerical integration of equation H5.20|l . witli tlie asymptotic 
result for the same values of N = 3, 10, 23, 40 and 1 < L < 100 as in FigureO 

The number of significant m = eigenvalues can be even more simply approximated 
by Nq w 2\/]V/7r w (L + l)Q/n, as shown. 




9 m 



w= to 




5 




(N=.23j 










I 1 

if ■ ■♦• 










.r- 













.9 2 4 B 8 10 




Fig. 7.2. Comparison of the nu mber N m of significant eigenvalues affixed order m (grey) 
with the asymptotic approximation {7.2^ (black). The Shannon number N — 3, 10, 23, 40 is 
kept constant in each of the four panels, and the bandwidth used to compute the exact values 
of Nm varies between L = 1 (worst fitting) and L = 100 (best fitting). Points inconsistent 
with the constraint A/{4'k) = N/{L + 1)^ < 1 are not plotted. White triangles show the 
simplified zonal (m = 0^ approximation No ~ (L + 1)0/tt. 



8. Conclusion. An orthogonal family of bandlimited spherical harmonic expan- 
sions that are optimally concentrated within a finite region R of the unit sphere can 
be computed by solving either a symmetric matrix eigenvalue problem in the spectral 
domain or an equivalent Fredholm integral eigenvalue problem in the spatial domain. 
Every eigenvalue < A < 1 is a measure of both the spatial concentration of the 
bandlimited eigenfunction g(f) and the spectral concentration of the spacelimited 
eigenfunction h{r) that coincides with g(r) inside the region of concentration. The 
number of well concentrated eigenfunctions is = (L + 1)'^A/ (An), where L is the 
maximal spherical degree and A is the area of the region of concentration. Roughly 
speaking, this Shannon number N is the dimension of the space of functions /(f) that 
can be concentrated both within a finite region R of the sphere and within a spectral 
interval < ^ < L. For a small region A <C 47r, and a moderate maximal spherical 
harmonic degree L, the optimally concentrated bandlimited eigenfunctions g{r) and 
associated spacelimited eigenfunctions h{r) can be computed accurately, even for an 
irregularly shaped region R. In the special, but important, case of a circular polar 
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cap, every eigcnfunction can be computed accurately, by numerical diagonalization 
of a commuting tridiagonal matrix, which has a simple Sturm-Liouvillc spectrum. 
Just as Slepian's onc-dimcnsional. prolate spheroidal eigentapcrs have proven to be 
extremely useful in time-frequency spectral analysis, we expect the two-dimensional, 
spherical eigenfunctions developed here to have a wide variety of spatiospectral data 
analysis applications in fields such as geophysics, planetary science and cosmology. 
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Appendix A. Computational considerations. Here we present a brief de- 
scription of the numerical methods employed in this study. All of the computations 
described here have been performed using double precision arithmetic. Statements 
regarding machine precision refer to double precision, with a roundoff error of ~ 10^^^. 

A.l. Concentration within a polar cap. We compute the colatitudinal eigen- 
functions gi{6), g2{0), ■ ■ ■ , gL-m+i{0) of an axisymmmetric polar cap < 6* < using 
three different methods. The first method is by numerical diagonalization of the 
(L — TO + 1) X (L — TO -1- 1) matrix D in equation (|5.4|l . We do not implement the 
Wigner 3-j expressio n (15.711 for the matrix elements Dw , but instead use Gauss- 
Legendre quadrature |43 to evaluate the defining integral H5.6|l : 

Dw = I A'i,„(arccos/x)X;/„i(arccos/x) 

cos 
,7 

(A.l) w y^WjXimi&rccos fij)Xi'r,t{arccosfij), 

where /ii,yU2, ■ ■ ■ ^f-J ^-re the roots of the Legendre polynomial Pj{fl), rescaled from 
-1 < Pj < 1 to cose < ^iJ < 1, and Wj = 2(1 - /if)"M-PJ(Mi)]"^ i = 1, 2, . . . , J are 
the associated integration weights. Only the uppermost triangular matrix elements 
Dill ,1^1' are computed explicitly; the lowermost elements are infilled using the 
symmetry Dui ~ Di'i. The order J of the Gauss-Legendre integration is adjusted 
upward mitil the i — to + 1 spatial-domain eigenfunctions gi (6) , g2 (0) , . . . , gL-m+i (^) 
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satisfy the orthogonality relations H5.11|l to within machine precision. The same high- 
order Gauss-Legendre quadrature rule is used to evaluate the orthogonality integrals. 
The Legcndre functions Xim (6) are computed with high accuracy to very high degree 
(I w 500) using a recursive algorithm 34iE3- 

The second method is by numerical solution of the Fredholm equation (|5.17b|) . 
Using a Gauss-Legendre quadrature rule to discretize this equation, we obtain 

J 

(A.2) ^u;,Z?(/.,,^;X/i;-) - A/i(/.,), J = 1,2,..., J. 

Equation HA.2|) can be rewritten as a symmetric algebraic eigenvalue equation, 
(A.3) (W^''^DW^/^)(W^/^h) = A(W^/^h), 

where h is a J-dimensional column vector with entries hj = h{fij), and where D 
and W denote the J x J matrices with elements Djji = D{^j,ij!^) and Wjji = 
WjSjji. The eigenvalues A and transformed eigenvectors W^'^h are computed by 
numerical diagonalization of the matrix W^'^DW^''^. The order of integration J is 
again chosen to ensure accurate orthogonality of the spatial-domain eigenfunctions 
hi{6), h2{0), ■ ■ ■ , hL^rn+i{0)- In the zonal (m = 0) case the choice J = L + 1 renders 
both of the integrations (|A.ip and (|A.2|I exact; for m Owe use a conservative, larger 
integration order J, since the integrands are no longer polynomials. The fixed-order 
power spectra shown in Figure |^| are computed by transforming the spatial-domain 
eigenfunctions of equation (|4.19|) to the spectral domain, using a spherical harmonic 
degree range m < I < 127 sufficient to avoid ahasing [HI ]. 

Even for moderate values of the maximal degree L and cap radius 8, the smallest 
eigenvalues . . . , XL-m, XL-m+i fall below machine precision. The associated, least well 
concentrated eigenfunctions computed using either of the above two direct methods 
are in that case essentially arbitrary orthogonal members of a numerically degenerate 
eigenspace, and are no longer accurate 0. Because of this, it is not possible to find 
the optimally excluded eigenfunctions of a small polar cap, or equivalently the opti- 
mally concentrated eigenfunctions of a large cap, by diagonalization of either of the 
matrices D or D. Fortunately, this difficulty can be overcome by the third method, 
which is numerical diagonalization of the tridiagonal Griinbaum matrix H5.44|) . The 
roughly equant spacing of the Sturm-Liouville eigenvalues Xi, X2, ■ • ■ , Xi-m+i enables 
all of the associated eigenfunctions to be calculated to within machine precision. The 
spatiospectral concentration factors Ai, A2, . . . , XL-m+i are computed to the same pre- 
cision, either by a posteriori matrix multiplication, A = g^Dg, or by Gauss-Legendre 
integration of the equivalent spatial relation (|5.11() . Both the significant and the in- 
significant eigenvalues computed using each of the above methods agree to within 
machine precision, providing a useful numerical check. Diagonalization of the tridiag- 
onal matrix G is the only numerically stable way to solve the concentration problem 
for either a large polar cap or a large maximal degree L. By an extension of the above 
analysis, it is even possible to use the Griinbaum operator Q to compute spacelimited 
eigenfunctions /iL-m+2(^), ^L-m+3(^), • ■ • that are in the null space 

A.2. Concentration within an arbitrarily shaped region. We solve the 
spatiospectral concentration problem for an arbitrarily shaped region R by numerical 
diagonalization of the {L + 1)^ x [L + 1)^ matrix D, with elements Dim,i'm' defined 
by equation H4.4|l . Given the splined boundary of R, we first find the northernmost 
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and southernmost points, with colatitudcs 0n and dg. For every 0^ < < Og, we 
then find the easternmost and westernmost points, with longitudes 4>c{()) and 4)^{9). 
In the case of a non-convex region with indentations and protuberances, there may 
be several such easternmost and westernmost points, which we shall index with an 
additional subscript i — 1,2, ... ,1. The integral over longitude, 



(A.4) ^m,n'{0)=J2 



cos m(j) 
sin m(p I 1 sinm' 



can be done analytically, and we use Gauss-Legendre quadrature to compute the 
remaining integral over colatitude: 

Dim,i'm' = / X;m(arccos/x)X;/„i/(arccos/x)$mm/(arccos/^) d/i 
J 

(A.5) « WjXimj&rCCOS flj)Xi>m' (arCCOS flj)^mm' (arCCOS flj). 

As in the case of a polar cap, we adjust the order of the integration J upward until 
the spatial-domain eigcnfunctions 5i(f ), (72(f),---, <7(l+i)2 (f ) satisfy the orthogonality 
relations (|4.11|) to within machine precision. There is no analogue of the Griinbaum 
operator Q for an arbitrarily shaped region, so only the eigcnfunctions associated with 
eigenvalues that are above machine precision can be computed accurately. In most 
practical applications, this is not a limitation, since we are generally interested only 
in the computable, well concentrated eigcnfunctions gi(f), 52(f), . . . ,gAr(f), which are 
associated with the numerically significant eigenvalues Ai, A2, . . . , A^v- 

A.3. Concentration within a non-polar circular cap. One of the principal 
applications of spherical Slepian functions in geophysics and planetary physics will 
be to analyze measurements within a circularl y sy mmetric region centered upon an 
arbitrary geographical location 9o, 0o [e-g-,|2a ISa; |51( |52l- The preferred procedure 
for determining the required optimally concentrated eigcnfunctions is first to com- 
pute the spherical harmonic coefficients gim of the eigcnfunctions H5.12|l concentrated 
within a pola r cap < 9 < 9, and then to rotate these to the desired cap location 
0; 0; fl3 I37II . The actual windowing of the data for further analysis may either be 
carried out in the spectral domain |52l |. or, more simply, by straightforward multipli- 
cation after transformation of the rotated eigcnfunctions to the spatial domain. If one 
wishes to avoid spherical harmonic rotation, it is also possible to compute the rotated 
eigcnfunctions directly, by performing the numerical integration in equation ljA.4|l on 
the analytically prescribed boundary of a cap of radius & centered at 9o,4>q, given by 

/ A „N , , A , A ,//,N arccos(cos9 — cos^cos^o) 
(A.6) (/)w,o(6') = </'o T A0(6') where A(t){9) = ^ 



sm u sm 



Appendix B. Spatiospectral concentration on a plane. In one of his many 
papers extending the one-dimensional analysis, Slepian js^ considered the spatiospec- 
tral concentration problem in a Cartesian space of arbitrary dimension. We present 
a brief review of the two-dimensional Cartesian concentration problem here, for com- 
parison with the flat-earth asymptotic analysis presented in Sectional 
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An arbitrary, rcal-valucd, squarc-intcgrablc function /(r) on the plane has the 
two-dimensional Fourier representation, analogous to the spherical harmonic repre- 
sentation 

/ -1 \ 2 /.oo POO 

(B.l) f{r)^i^—jj F(k)e*'^ "-d2k, H^) = J /We-^"^ "- d^r. 
The Fourier orthonormality relation analogous to equation H3.6|l is 

(B.2) ( ^ V r e^'^-f-^-'-') d'k = 6(r r') = ^^1^—^. 

V2^/i-oo 2^||r-r'|| 

Parseval's relation stipulates that the power of any function /(r) in the spatial and 
spectral domains is identical: 

Equation ljB.3|) is the planar analogue of the spherical relation ||/||^ = |lf|l^- 
We use g(r) to denote a bandlimited function, 



(B.4) 5(r)=(^) / G(k)e-"-d^k, 



1 



\M<K 



with no power above a maximal wavenumber K. By analogy with the optimization 
criterion 14.2|l . we seek to concentrate the power of g{r) into a finite region R: 

g d r 

(B.5) A = = maximum. 

9'd\ 



Bandlimited functions g{r) that maximize the ratio A in equation ljB.5p arc solutions 
to the Fourier domain eigenvalue equation, analogous to equation l|4.7|l . 

(B.6) f £»(k,k')G(k)(i2k' = AG(k), ||k|| < if, 

•^||k||<if 

where 

(B.7) D{-k,-k') ^ Q-^ J e''^^-^'^'' d^r. 

The corresponding problem in the spatial domain, analogous to equation H4.15|l . is 
(B.8) / i?(r,r')5(r)rfV = Ag(r), |r| < oo, 



where 



(B.9) D{r,r')=(^^f e'^ '^"-''"' d^k. 
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Spacclimited eigenfunctions h{r), which vanish outside of the region R, satisfy the 
same eigenvalue equation ljB.8|) . but with the domain of solution properly restricted: 



(B.IO) [ D{r,r')h{r)d^r' = \h{r), r e R. 

Jr 

The associated eigenvalue < A < 1 is a measure of both the spatial concentration of 
g(r) within the region R and the spectral concentration of h{r) within the wavcnumbcr 
range ||k|| < K. 

For consistency with H4.9|l . we rank order the eigenvalues so that Ai > A2 > . . .. 
The bandlimited, spatial-domain eigenfunctions gi(r), (7i(r), . . . may be chosen to be 
orthonormal over the whole plane ||r|| < 00 and orthogonal over the region R: 

(B.ll) / gaQpd^r ^ 6af3 and / ga9f3 <Pr = K^ap- 

J -00 J R 

Ths sum of the eigenvalues, or Shannon number, is given by 

°° r r A 

(B.12) iV = ^A„= / D{k,'k)(fk^ / D{r,r)d^r ^ —, 



where A is the area of the region of concentration R. Equations (|B.11|) and (|B.12|I 
are the planar analogues of the spherical relations H4.11|l and (|4.28|) . 

Comparison of equations (|B.2|) and IjB.Qp shows that the spatial-domain kernel 
D(r,r'), like its spherical counterpart D(f,f'), is a bandlimited delta function. Upon 
introducing polar coordinates and integrating over the angle, we can reduce D{r,r') 
to a form reminiscent of the representation H7.7|) : 

KMK\\r-r'\\) 



(B.13) D{v,v')^^[ Mk\\v-v'\\)kdk 



27r||r-r'|| 

Upon introducing scaled independent and dependent variables analogous to (|7.9|l . 

/ 47r / 47r 

(B.14) x=-^— r, x' = ^— r', .^(x) - .g(r), ^/.(x') - .g(r'), 

we can rewrite equations (|B.8|I and (|B.13p in a form identical to (|7.11|) - (|7.12|l and 
analogous to (|2.6() : 

,B.5) « f ..(A1|x .11) 

27r./B x-x 



where i?* is the image of the region of concentration R under the mapping (jB.14|l . 

If the region of concentration i? is a circle of radius Q, then a polar coordinate, 
r =: (g, (/)), representation analogous to (|5.12|) . 

\/2 g{q) coarrKj) if — L < to < 
(B.16) g{q,4>) = { g{q) if to = 

\/2 g{q) smmcj) if < 771 < L, 
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may be used to decompose equations (IB.8|) and (jB.13|) into a series of fixed-order 
eigenvalue problems analogous to H7.f 7|I - H7.f 8|l : 



(B.f7) / D{q,q')g{q')q'dq' ^\g{q), 

Jo 

where 

(B.f8) D{q,q')=K^[ J^{Kpq)J^{Kpq')pdp. 

Jo 

The transformations 

(B.f9) x = (z/g, x'=q'/Q, ^{x)=g{q), i^{x') = g{q') 



convert equations ljB.17|) - (|B.18p into the scaled eigenvalue problem 
(B.20) 4Af /" / J,n{2yfNpx).ln{2yfNpx')pdpi:{x')x'dx' = Xi){x), 



Jo 



which is identical to (I7j0|) - (|7.21() . and dependent only upon the Shannon number 
N = jK^Q^. Slepian |52| has noted that equation ljB.20|l is an iterated version of the 
equivalent "square root" equation 



(B.2f) 2Vn J,n{2y/Nxx')i:{x')x'dx' = V\^ix) 

Jo 



The eigenvalues Ai, A2, . . . and eigenfunctions "01(2;), '02 (a^), ■ • • of equation ljB.20|) may 
alternatively be found by solving the equivalent equation l|B.21|) . 

In the asymptotic limit (|7.1|) . both the general and axisymmetric spherical con- 
centration problems are seen to be identical to the corresponding concentration prob- 
lem in a plane, with the maximal wavenumbcr K replaced by the integer L + I. 
The planar problem exhibits exact Shannon-number scaling analogous to that of the 
one-dimensional problem H2.6|l . whereas the scaling of the spherical problem is only 
asymptotic. Equation H7.22(l giving the number of significant eigenvalues Nm associ- 
ated with each angular order m is exact in the case of the plane. 



